home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  26.8 KB  |  931 lines
  1. /* specfunc/bessel.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21. /* Miscellaneous support functions for Bessel function evaluations.
  22.  */
  23. #include <config.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_errno.h>
  26. #include <gsl/gsl_sf_airy.h>
  27. #include <gsl/gsl_sf_elementary.h>
  28. #include <gsl/gsl_sf_exp.h>
  29. #include <gsl/gsl_sf_gamma.h>
  30. #include <gsl/gsl_sf_trig.h>
  31.  
  32. #include "error.h"
  33.  
  34. #include "bessel_amp_phase.h"
  35. #include "bessel_temme.h"
  36. #include "bessel.h"
  37.  
  38. #define CubeRoot2_  1.25992104989487316476721060728
  39.  
  40.  
  41.  
  42. /* Debye functions [Abramowitz+Stegun, 9.3.9-10] */
  43.  
  44. inline static double 
  45. debye_u1(const double * tpow)
  46. {
  47.   return (3.0*tpow[1] - 5.0*tpow[3])/24.0;
  48. }
  49.  
  50. inline static double 
  51. debye_u2(const double * tpow)
  52. {
  53.   return (81.0*tpow[2] - 462.0*tpow[4] + 385.0*tpow[6])/1152.0;
  54. }
  55.  
  56. inline
  57. static double debye_u3(const double * tpow)
  58. {
  59.   return (30375.0*tpow[3] - 369603.0*tpow[5] + 765765.0*tpow[7] - 425425.0*tpow[9])/414720.0;
  60. }
  61.  
  62. inline
  63. static double debye_u4(const double * tpow)
  64. {
  65.   return (4465125.0*tpow[4] - 94121676.0*tpow[6] + 349922430.0*tpow[8] - 
  66.           446185740.0*tpow[10] + 185910725.0*tpow[12])/39813120.0;
  67. }
  68.  
  69. inline
  70. static double debye_u5(const double * tpow)
  71. {
  72.   return (1519035525.0*tpow[5]     - 49286948607.0*tpow[7] + 
  73.           284499769554.0*tpow[9]   - 614135872350.0*tpow[11] + 
  74.           566098157625.0*tpow[13]  - 188699385875.0*tpow[15])/6688604160.0;
  75. }
  76.  
  77. #if 0
  78. inline
  79. static double debye_u6(const double * tpow)
  80. {
  81.   return (2757049477875.0*tpow[6] - 127577298354750.0*tpow[8] + 
  82.           1050760774457901.0*tpow[10] - 3369032068261860.0*tpow[12] + 
  83.           5104696716244125.0*tpow[14] - 3685299006138750.0*tpow[16] + 
  84.           1023694168371875.0*tpow[18])/4815794995200.0;
  85. }
  86. #endif
  87.  
  88.  
  89. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  90.  
  91. int
  92. gsl_sf_bessel_IJ_taylor_e(const double nu, const double x,
  93.                              const int sign,
  94.                              const int kmax,
  95.                              const double threshold,
  96.                              gsl_sf_result * result
  97.                              )
  98. {
  99.   /* CHECK_POINTER(result) */
  100.  
  101.   if(nu < 0.0 || x < 0.0) {
  102.     DOMAIN_ERROR(result);
  103.   }
  104.   else if(x == 0.0) {
  105.     if(nu == 0.0) {
  106.       result->val = 1.0;
  107.       result->err = 0.0;
  108.     }
  109.     else {
  110.       result->val = 0.0;
  111.       result->err = 0.0;
  112.     }
  113.     return GSL_SUCCESS;
  114.   }
  115.   else {
  116.     gsl_sf_result prefactor;   /* (x/2)^nu / Gamma(nu+1) */
  117.     gsl_sf_result sum;
  118.  
  119.     int stat_pre;
  120.     int stat_sum;
  121.     int stat_mul;
  122.  
  123.     if(nu == 0.0) {
  124.       prefactor.val = 1.0;
  125.       prefactor.err = 0.0;
  126.       stat_pre = GSL_SUCCESS;
  127.     }
  128.     else if(nu < INT_MAX-1) {
  129.       /* Separate the integer part and use
  130.        * y^nu / Gamma(nu+1) = y^N /N! y^f / (N+1)_f,
  131.        * to control the error.
  132.        */
  133.       const int    N = (int)floor(nu + 0.5);
  134.       const double f = nu - N;
  135.       gsl_sf_result poch_factor;
  136.       gsl_sf_result tc_factor;
  137.       const int stat_poch = gsl_sf_poch_e(N+1.0, f, &poch_factor);
  138.       const int stat_tc   = gsl_sf_taylorcoeff_e(N, 0.5*x, &tc_factor);
  139.       const double p = pow(0.5*x,f);
  140.       prefactor.val  = tc_factor.val * p / poch_factor.val;
  141.       prefactor.err  = tc_factor.err * p / poch_factor.val;
  142.       prefactor.err += fabs(prefactor.val) / poch_factor.val * poch_factor.err;
  143.       prefactor.err += 2.0 * GSL_DBL_EPSILON * fabs(prefactor.val);
  144.       stat_pre = GSL_ERROR_SELECT_2(stat_tc, stat_poch);
  145.     }
  146.     else {
  147.       gsl_sf_result lg;
  148.       const int stat_lg = gsl_sf_lngamma_e(nu+1.0, &lg);
  149.       const double term1  = nu*log(0.5*x);
  150.       const double term2  = lg.val;
  151.       const double ln_pre = term1 - term2;
  152.       const double ln_pre_err = GSL_DBL_EPSILON * (fabs(term1)+fabs(term2)) + lg.err;
  153.       const int stat_ex = gsl_sf_exp_err_e(ln_pre, ln_pre_err, &prefactor);
  154.       stat_pre = GSL_ERROR_SELECT_2(stat_ex, stat_lg);
  155.     }
  156.  
  157.     /* Evaluate the sum.
  158.      * [Abramowitz+Stegun, 9.1.10]
  159.      * [Abramowitz+Stegun, 9.6.7]
  160.      */
  161.     {
  162.       const double y = sign * 0.25 * x*x;
  163.       double sumk = 1.0;
  164.       double term = 1.0;
  165.       int k;
  166.  
  167.       for(k=1; k<=kmax; k++) {
  168.         term *= y/((nu+k)*k);
  169.         sumk += term;
  170.         if(fabs(term/sumk) < threshold) break;
  171.       }
  172.  
  173.       sum.val = sumk;
  174.       sum.err = threshold * fabs(sumk);
  175.  
  176.       stat_sum = ( k >= kmax ? GSL_EMAXITER : GSL_SUCCESS );
  177.     }
  178.  
  179.     stat_mul = gsl_sf_multiply_err_e(prefactor.val, prefactor.err,
  180.                                         sum.val, sum.err,
  181.                                         result);
  182.  
  183.     return GSL_ERROR_SELECT_3(stat_mul, stat_pre, stat_sum);
  184.   }
  185. }
  186.  
  187.  
  188. /* x >> nu*nu+1
  189.  * error ~ O( ((nu*nu+1)/x)^3 )
  190.  *
  191.  * empirical error analysis:
  192.  *   choose  GSL_ROOT3_MACH_EPS * x > (nu*nu + 1)
  193.  *
  194.  * This is not especially useful. When the argument gets
  195.  * large enough for this to apply, the cos() and sin()
  196.  * start loosing digits. However, this seems inevitable
  197.  * for this particular method.
  198.  */
  199. int
  200. gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result)
  201. {
  202.   double mu   = 4.0*nu*nu;
  203.   double mum1 = mu-1.0;
  204.   double mum9 = mu-9.0;
  205.   double chi = x - (0.5*nu + 0.25)*M_PI;
  206.   double P   = 1.0 - mum1*mum9/(128.0*x*x);
  207.   double Q   = mum1/(8.0*x);
  208.   double pre = sqrt(2.0/(M_PI*x));
  209.   double c   = cos(chi);
  210.   double s   = sin(chi);
  211.   double r   = mu/x;
  212.   result->val  = pre * (c*P - s*Q);
  213.   result->err  = pre * GSL_DBL_EPSILON * (fabs(c*P) + fabs(s*Q));
  214.   result->err += pre * fabs(0.1*r*r*r);
  215.   return GSL_SUCCESS;
  216. }
  217.  
  218.  
  219. /* x >> nu*nu+1
  220.  */
  221. int
  222. gsl_sf_bessel_Ynu_asympx_e(const double nu, const double x, gsl_sf_result * result)
  223. {
  224.   double ampl;
  225.   double theta;
  226.   double alpha = x;
  227.   double beta  = -0.5*nu*M_PI;
  228.   int stat_a = gsl_sf_bessel_asymp_Mnu_e(nu, x, &l);
  229.   int stat_t = gsl_sf_bessel_asymp_thetanu_corr_e(nu, x, &theta);
  230.   double sin_alpha = sin(alpha);
  231.   double cos_alpha = cos(alpha);
  232.   double sin_chi   = sin(beta + theta);
  233.   double cos_chi   = cos(beta + theta);
  234.   double sin_term     = sin_alpha * cos_chi + sin_chi * cos_alpha;
  235.   double sin_term_mag = fabs(sin_alpha * cos_chi) + fabs(sin_chi * cos_alpha);
  236.   result->val  = ampl * sin_term;
  237.   result->err  = fabs(ampl) * GSL_DBL_EPSILON * sin_term_mag;
  238.   result->err += fabs(result->val) * 2.0 * GSL_DBL_EPSILON;
  239.  
  240.   if(fabs(alpha) > 1.0/GSL_DBL_EPSILON) {
  241.     result->err *= 0.5 * fabs(alpha);
  242.   }
  243.   else if(fabs(alpha) > 1.0/GSL_SQRT_DBL_EPSILON) {
  244.     result->err *= 256.0 * fabs(alpha) * GSL_SQRT_DBL_EPSILON;
  245.   }
  246.  
  247.   return GSL_ERROR_SELECT_2(stat_t, stat_a);
  248. }
  249.  
  250.  
  251. /* x >> nu*nu+1
  252.  */
  253. int
  254. gsl_sf_bessel_Inu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
  255. {
  256.   double mu   = 4.0*nu*nu;
  257.   double mum1 = mu-1.0;
  258.   double mum9 = mu-9.0;
  259.   double pre  = 1.0/sqrt(2.0*M_PI*x);
  260.   double r    = mu/x;
  261.   result->val = pre * (1.0 - mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
  262.   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
  263.   return GSL_SUCCESS;
  264. }
  265.  
  266. /* x >> nu*nu+1
  267.  */
  268. int
  269. gsl_sf_bessel_Knu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
  270. {
  271.   double mu   = 4.0*nu*nu;
  272.   double mum1 = mu-1.0;
  273.   double mum9 = mu-9.0;
  274.   double pre  = sqrt(M_PI/(2.0*x));
  275.   double r    = nu/x;
  276.   result->val = pre * (1.0 + mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
  277.   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
  278.   return GSL_SUCCESS;
  279. }
  280.  
  281.  
  282. /* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.7]
  283.  *
  284.  * error:
  285.  *   The error has the form u_N(t)/nu^N  where  0 <= t <= 1.
  286.  *   It is not hard to show that |u_N(t)| is small for such t.
  287.  *   We have N=6 here, and |u_6(t)| < 0.025, so the error is clearly
  288.  *   bounded by 0.025/nu^6. This gives the asymptotic bound on nu
  289.  *   seen below as nu ~ 100. For general MACH_EPS it will be 
  290.  *                     nu > 0.5 / MACH_EPS^(1/6)
  291.  *   When t is small, the bound is even better because |u_N(t)| vanishes
  292.  *   as t->0. In fact u_N(t) ~ C t^N as t->0, with C ~= 0.1.
  293.  *   We write
  294.  *                     err_N <= min(0.025, C(1/(1+(x/nu)^2))^3) / nu^6
  295.  *   therefore
  296.  *                     min(0.29/nu^2, 0.5/(nu^2+x^2)) < MACH_EPS^{1/3}
  297.  *   and this is the general form.
  298.  *
  299.  * empirical error analysis, assuming 14 digit requirement:
  300.  *   choose   x > 50.000 nu   ==>  nu >   3
  301.  *   choose   x > 10.000 nu   ==>  nu >  15
  302.  *   choose   x >  2.000 nu   ==>  nu >  50
  303.  *   choose   x >  1.000 nu   ==>  nu >  75
  304.  *   choose   x >  0.500 nu   ==>  nu >  80
  305.  *   choose   x >  0.100 nu   ==>  nu >  83
  306.  *
  307.  * This makes sense. For x << nu, the error will be of the form u_N(1)/nu^N,
  308.  * since the polynomial term will be evaluated near t=1, so the bound
  309.  * on nu will become constant for small x. Furthermore, increasing x with
  310.  * nu fixed will decrease the error.
  311.  */
  312. int
  313. gsl_sf_bessel_Inu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
  314. {
  315.   int i;
  316.   double z = x/nu;
  317.   double root_term = sqrt(1.0 + z*z);
  318.   double pre = 1.0/sqrt(2.0*M_PI*nu * root_term);
  319.   double eta = root_term + log(z/(1.0+root_term));
  320.   double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(-z + eta) : -0.5*nu/z*(1.0 - 1.0/(12.0*z*z)) );
  321.   gsl_sf_result ex_result;
  322.   int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
  323.   if(stat_ex == GSL_SUCCESS) {
  324.     double t = 1.0/root_term;
  325.     double sum;
  326.     double tpow[16];
  327.     tpow[0] = 1.0;
  328.     for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
  329.     sum = 1.0 + debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) + debye_u3(tpow)/(nu*nu*nu)
  330.           + debye_u4(tpow)/(nu*nu*nu*nu) + debye_u5(tpow)/(nu*nu*nu*nu*nu);
  331.     result->val  = pre * ex_result.val * sum;
  332.     result->err  = pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
  333.     result->err += pre * ex_result.err * fabs(sum);
  334.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  335.     return GSL_SUCCESS;
  336.   }
  337.   else {
  338.     result->val = 0.0;
  339.     result->err = 0.0;
  340.     return stat_ex;
  341.   }
  342. }
  343.  
  344.  
  345. /* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.8]
  346.  *
  347.  * error:
  348.  *   identical to that above for Inu_scaled
  349.  */
  350. int
  351. gsl_sf_bessel_Knu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
  352. {
  353.   int i;
  354.   double z = x/nu;
  355.   double root_term = sqrt(1.0 + z*z);
  356.   double pre = sqrt(M_PI/(2.0*nu*root_term));
  357.   double eta = root_term + log(z/(1.0+root_term));
  358.   double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(z - eta) : 0.5*nu/z*(1.0 + 1.0/(12.0*z*z)) );
  359.   gsl_sf_result ex_result;
  360.   int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
  361.   if(stat_ex == GSL_SUCCESS) {
  362.     double t = 1.0/root_term;
  363.     double sum;
  364.     double tpow[16];
  365.     tpow[0] = 1.0;
  366.     for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
  367.     sum = 1.0 - debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) - debye_u3(tpow)/(nu*nu*nu)
  368.           + debye_u4(tpow)/(nu*nu*nu*nu) - debye_u5(tpow)/(nu*nu*nu*nu*nu);
  369.     result->val  = pre * ex_result.val * sum;
  370.     result->err  = pre * ex_result.err * fabs(sum);
  371.     result->err += pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
  372.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  373.     return GSL_SUCCESS;
  374.   }
  375.   else {
  376.     result->val = 0.0;
  377.     result->err = 0.0;
  378.     return stat_ex;
  379.   }
  380. }
  381.  
  382.  
  383. /* Evaluate J_mu(x),J_{mu+1}(x) and Y_mu(x),Y_{mu+1}(x)  for |mu| < 1/2
  384.  */
  385. int
  386. gsl_sf_bessel_JY_mu_restricted(const double mu, const double x,
  387.                                gsl_sf_result * Jmu, gsl_sf_result * Jmup1,
  388.                                gsl_sf_result * Ymu, gsl_sf_result * Ymup1)
  389. {
  390.   /* CHECK_POINTER(Jmu) */
  391.   /* CHECK_POINTER(Jmup1) */
  392.   /* CHECK_POINTER(Ymu) */
  393.   /* CHECK_POINTER(Ymup1) */
  394.  
  395.   if(x < 0.0 || fabs(mu) > 0.5) {
  396.     Jmu->val   = 0.0;
  397.     Jmu->err   = 0.0;
  398.     Jmup1->val = 0.0;
  399.     Jmup1->err = 0.0;
  400.     Ymu->val   = 0.0;
  401.     Ymu->err   = 0.0;
  402.     Ymup1->val = 0.0;
  403.     Ymup1->err = 0.0;
  404.     GSL_ERROR ("error", GSL_EDOM);
  405.   }
  406.   else if(x == 0.0) {
  407.     if(mu == 0.0) {
  408.       Jmu->val   = 1.0;
  409.       Jmu->err   = 0.0;
  410.     }
  411.     else {
  412.       Jmu->val   = 0.0;
  413.       Jmu->err   = 0.0;
  414.     }
  415.     Jmup1->val = 0.0;
  416.     Jmup1->err = 0.0;
  417.     Ymu->val   = 0.0;
  418.     Ymu->err   = 0.0;
  419.     Ymup1->val = 0.0;
  420.     Ymup1->err = 0.0;
  421.     GSL_ERROR ("error", GSL_EDOM);
  422.   }
  423.   else {
  424.     int stat_Y;
  425.     int stat_J;
  426.  
  427.     if(x < 2.0) {
  428.       /* Use Taylor series for J and the Temme series for Y.
  429.        * The Taylor series for J requires nu > 0, so we shift
  430.        * up one and use the recursion relation to get Jmu, in
  431.        * case mu < 0.
  432.        */
  433.       gsl_sf_result Jmup2;
  434.       int stat_J1 = gsl_sf_bessel_IJ_taylor_e(mu+1.0, x, -1, 100, GSL_DBL_EPSILON,  Jmup1);
  435.       int stat_J2 = gsl_sf_bessel_IJ_taylor_e(mu+2.0, x, -1, 100, GSL_DBL_EPSILON, &Jmup2);
  436.       double c = 2.0*(mu+1.0)/x;
  437.       Jmu->val  = c * Jmup1->val - Jmup2.val;
  438.       Jmu->err  = c * Jmup1->err + Jmup2.err;
  439.       Jmu->err += 2.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
  440.       stat_J = GSL_ERROR_SELECT_2(stat_J1, stat_J2);
  441.       stat_Y = gsl_sf_bessel_Y_temme(mu, x, Ymu, Ymup1);
  442.       return GSL_ERROR_SELECT_2(stat_J, stat_Y);
  443.     }
  444.     else if(x < 1000.0) {
  445.       double P, Q;
  446.       double J_ratio;
  447.       double J_sgn;
  448.       const int stat_CF1 = gsl_sf_bessel_J_CF1(mu, x, &J_ratio, &J_sgn);
  449.       const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
  450.       double Jprime_J_ratio = mu/x - J_ratio;
  451.       double gamma = (P - Jprime_J_ratio)/Q;
  452.       Jmu->val = J_sgn * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jprime_J_ratio)));
  453.       Jmu->err = 4.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
  454.       Jmup1->val = J_ratio * Jmu->val;
  455.       Jmup1->err = fabs(J_ratio) * Jmu->err;
  456.       Ymu->val = gamma * Jmu->val;
  457.       Ymu->err = fabs(gamma) * Jmu->err;
  458.       Ymup1->val = Ymu->val * (mu/x - P - Q/gamma);
  459.       Ymup1->err = Ymu->err * fabs(mu/x - P - Q/gamma) + 4.0*GSL_DBL_EPSILON*fabs(Ymup1->val);
  460.       return GSL_ERROR_SELECT_2(stat_CF1, stat_CF2);
  461.     }
  462.     else {
  463.       /* Use asymptotics for large argument.
  464.        */
  465.       const int stat_J0 = gsl_sf_bessel_Jnu_asympx_e(mu,     x, Jmu);
  466.       const int stat_J1 = gsl_sf_bessel_Jnu_asympx_e(mu+1.0, x, Jmup1);
  467.       const int stat_Y0 = gsl_sf_bessel_Ynu_asympx_e(mu,     x, Ymu);
  468.       const int stat_Y1 = gsl_sf_bessel_Ynu_asympx_e(mu+1.0, x, Ymup1);
  469.       stat_J = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
  470.       stat_Y = GSL_ERROR_SELECT_2(stat_Y0, stat_Y1);
  471.       return GSL_ERROR_SELECT_2(stat_J, stat_Y);
  472.     }
  473.   }
  474. }
  475.  
  476.  
  477. int
  478. gsl_sf_bessel_J_CF1(const double nu, const double x,
  479.                     double * ratio, double * sgn)
  480. {
  481.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  482.   const int maxiter = 10000;
  483.   int n = 1;
  484.   double Anm2 = 1.0;
  485.   double Bnm2 = 0.0;
  486.   double Anm1 = 0.0;
  487.   double Bnm1 = 1.0;
  488.   double a1 = x/(2.0*(nu+1.0));
  489.   double An = Anm1 + a1*Anm2;
  490.   double Bn = Bnm1 + a1*Bnm2;
  491.   double an;
  492.   double fn = An/Bn;
  493.   double dn = a1;
  494.   double s  = 1.0;
  495.  
  496.   while(n < maxiter) {
  497.     double old_fn;
  498.     double del;
  499.     n++;
  500.     Anm2 = Anm1;
  501.     Bnm2 = Bnm1;
  502.     Anm1 = An;
  503.     Bnm1 = Bn;
  504.     an = -x*x/(4.0*(nu+n-1.0)*(nu+n));
  505.     An = Anm1 + an*Anm2;
  506.     Bn = Bnm1 + an*Bnm2;
  507.  
  508.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  509.       An /= RECUR_BIG;
  510.       Bn /= RECUR_BIG;
  511.       Anm1 /= RECUR_BIG;
  512.       Bnm1 /= RECUR_BIG;
  513.       Anm2 /= RECUR_BIG;
  514.       Bnm2 /= RECUR_BIG;
  515.     }
  516.  
  517.     old_fn = fn;
  518.     fn = An/Bn;
  519.     del = old_fn/fn;
  520.  
  521.     dn = 1.0 / (2.0*(nu+n)/x - dn);
  522.     if(dn < 0.0) s = -s;
  523.  
  524.     if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
  525.   }
  526.  
  527.   *ratio = fn;
  528.   *sgn   = s;
  529.  
  530.   if(n >= maxiter)
  531.     GSL_ERROR ("error", GSL_EMAXITER);
  532.   else
  533.     return GSL_SUCCESS;
  534. }
  535.  
  536.  
  537.  
  538. /* Evaluate the continued fraction CF1 for J_{nu+1}/J_nu
  539.  * using Gautschi (Euler) equivalent series.
  540.  * This exhibits an annoying problem because the
  541.  * a_k are not positive definite (in fact they are all negative).
  542.  * There are cases when rho_k blows up. Example: nu=1,x=4.
  543.  */
  544. #if 0
  545. int
  546. gsl_sf_bessel_J_CF1_ser(const double nu, const double x,
  547.                         double * ratio, double * sgn)
  548. {
  549.   const int maxk = 20000;
  550.   double tk   = 1.0;
  551.   double sum  = 1.0;
  552.   double rhok = 0.0;
  553.   double dk = 0.0;
  554.   double s  = 1.0;
  555.   int k;
  556.  
  557.   for(k=1; k<maxk; k++) {
  558.     double ak = -0.25 * (x/(nu+k)) * x/(nu+k+1.0);
  559.     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
  560.     tk  *= rhok;
  561.     sum += tk;
  562.  
  563.     dk = 1.0 / (2.0/x - (nu+k-1.0)/(nu+k) * dk);
  564.     if(dk < 0.0) s = -s;
  565.  
  566.     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  567.   }
  568.  
  569.   *ratio = x/(2.0*(nu+1.0)) * sum;
  570.   *sgn   = s;
  571.  
  572.   if(k == maxk)
  573.     GSL_ERROR ("error", GSL_EMAXITER);
  574.   else
  575.     return GSL_SUCCESS;
  576. }
  577. #endif
  578.  
  579.  
  580. /* Evaluate the continued fraction CF1 for I_{nu+1}/I_nu
  581.  * using Gautschi (Euler) equivalent series.
  582.  */
  583. int
  584. gsl_sf_bessel_I_CF1_ser(const double nu, const double x, double * ratio)
  585. {
  586.   const int maxk = 20000;
  587.   double tk   = 1.0;
  588.   double sum  = 1.0;
  589.   double rhok = 0.0;
  590.   int k;
  591.  
  592.   for(k=1; k<maxk; k++) {
  593.     double ak = 0.25 * (x/(nu+k)) * x/(nu+k+1.0);
  594.     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
  595.     tk  *= rhok;
  596.     sum += tk;
  597.     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  598.   }
  599.  
  600.   *ratio = x/(2.0*(nu+1.0)) * sum;
  601.  
  602.   if(k == maxk)
  603.     GSL_ERROR ("error", GSL_EMAXITER);
  604.   else
  605.     return GSL_SUCCESS;
  606. }
  607.  
  608.  
  609. int
  610. gsl_sf_bessel_JY_steed_CF2(const double nu, const double x,
  611.                            double * P, double * Q)
  612. {
  613.   const int max_iter = 10000;
  614.   const double SMALL = 1.0e-100;
  615.  
  616.   int i = 1;
  617.  
  618.   double x_inv = 1.0/x;
  619.   double a = 0.25 - nu*nu;
  620.   double p = -0.5*x_inv;
  621.   double q = 1.0;
  622.   double br = 2.0*x;
  623.   double bi = 2.0;
  624.   double fact = a*x_inv/(p*p + q*q);
  625.   double cr = br + q*fact;
  626.   double ci = bi + p*fact;
  627.   double den = br*br + bi*bi;
  628.   double dr = br/den;
  629.   double di = -bi/den;
  630.   double dlr = cr*dr - ci*di;
  631.   double dli = cr*di + ci*dr;
  632.   double temp = p*dlr - q*dli;
  633.   q = p*dli + q*dlr;
  634.   p = temp;
  635.   for (i=2; i<=max_iter; i++) {
  636.     a  += 2*(i-1);
  637.     bi += 2.0;
  638.     dr = a*dr + br;
  639.     di = a*di + bi;
  640.     if(fabs(dr)+fabs(di) < SMALL) dr = SMALL;
  641.     fact = a/(cr*cr+ci*ci);
  642.     cr = br + cr*fact;
  643.     ci = bi - ci*fact;
  644.     if(fabs(cr)+fabs(ci) < SMALL) cr = SMALL;
  645.     den = dr*dr + di*di;
  646.     dr /= den;
  647.     di /= -den;
  648.     dlr = cr*dr - ci*di;
  649.     dli = cr*di + ci*dr;
  650.     temp = p*dlr - q*dli;
  651.     q = p*dli + q*dlr;
  652.     p = temp;
  653.     if(fabs(dlr-1.0)+fabs(dli) < GSL_DBL_EPSILON) break;
  654.   }
  655.  
  656.   *P = p;
  657.   *Q = q;
  658.  
  659.   if(i == max_iter)
  660.     GSL_ERROR ("error", GSL_EMAXITER);
  661.   else
  662.     return GSL_SUCCESS;
  663. }
  664.  
  665.  
  666. /* Evaluate continued fraction CF2, using Thompson-Barnett-Temme method,
  667.  * to obtain values of exp(x)*K_nu and exp(x)*K_{nu+1}.
  668.  *
  669.  * This is unstable for small x; x > 2 is a good cutoff.
  670.  * Also requires |nu| < 1/2.
  671.  */
  672. int
  673. gsl_sf_bessel_K_scaled_steed_temme_CF2(const double nu, const double x,
  674.                                        double * K_nu, double * K_nup1,
  675.                                        double * Kp_nu)
  676. {
  677.   const int maxiter = 10000;
  678.  
  679.   int i = 1;
  680.   double bi = 2.0*(1.0 + x);
  681.   double di = 1.0/bi;
  682.   double delhi = di;
  683.   double hi    = di;
  684.  
  685.   double qi   = 0.0;
  686.   double qip1 = 1.0;
  687.  
  688.   double ai = -(0.25 - nu*nu);
  689.   double a1 = ai;
  690.   double ci = -ai;
  691.   double Qi = -ai;
  692.  
  693.   double s = 1.0 + Qi*delhi;
  694.  
  695.   for(i=2; i<=maxiter; i++) {
  696.     double dels;
  697.     double tmp;
  698.     ai -= 2.0*(i-1);
  699.     ci  = -ai*ci/i;
  700.     tmp  = (qi - bi*qip1)/ai;
  701.     qi   = qip1;
  702.     qip1 = tmp;
  703.     Qi += ci*qip1;
  704.     bi += 2.0;
  705.     di  = 1.0/(bi + ai*di);
  706.     delhi = (bi*di - 1.0) * delhi;
  707.     hi += delhi;
  708.     dels = Qi*delhi;
  709.     s += dels;
  710.     if(fabs(dels/s) < GSL_DBL_EPSILON) break;
  711.   }
  712.   
  713.   hi *= -a1;
  714.   
  715.   *K_nu   = sqrt(M_PI/(2.0*x)) / s;
  716.   *K_nup1 = *K_nu * (nu + x + 0.5 - hi)/x;
  717.   *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
  718.   if(i == maxiter)
  719.     GSL_ERROR ("error", GSL_EMAXITER);
  720.   else
  721.     return GSL_SUCCESS;
  722. }
  723.  
  724.  
  725. int gsl_sf_bessel_cos_pi4_e(double y, double eps, gsl_sf_result * result)
  726. {
  727.   const double sy = sin(y);
  728.   const double cy = cos(y);
  729.   const double s = sy + cy;
  730.   const double d = sy - cy;
  731.   const double abs_sum = fabs(cy) + fabs(sy);
  732.   double seps;
  733.   double ceps;
  734.   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
  735.     const double e2 = eps*eps;
  736.     seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
  737.     ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
  738.   }
  739.   else {
  740.     seps = sin(eps);
  741.     ceps = cos(eps);
  742.   }
  743.   result->val = (ceps * s - seps * d)/ M_SQRT2;
  744.   result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;
  745.  
  746.   /* Try to account for error in evaluation of sin(y), cos(y).
  747.    * This is a little sticky because we don't really know
  748.    * how the library routines are doing their argument reduction.
  749.    * However, we will make a reasonable guess.
  750.    * FIXME ?
  751.    */
  752.   if(y > 1.0/GSL_DBL_EPSILON) {
  753.     result->err *= 0.5 * y;
  754.   }
  755.   else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
  756.     result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
  757.   }
  758.  
  759.   return GSL_SUCCESS;
  760. }
  761.  
  762.  
  763. int gsl_sf_bessel_sin_pi4_e(double y, double eps, gsl_sf_result * result)
  764. {
  765.   const double sy = sin(y);
  766.   const double cy = cos(y);
  767.   const double s = sy + cy;
  768.   const double d = sy - cy;
  769.   const double abs_sum = fabs(cy) + fabs(sy);
  770.   double seps;
  771.   double ceps;
  772.   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
  773.     const double e2 = eps*eps;
  774.     seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
  775.     ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
  776.   }
  777.   else {
  778.     seps = sin(eps);
  779.     ceps = cos(eps);
  780.   }
  781.   result->val = (ceps * d + seps * s)/ M_SQRT2;
  782.   result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;
  783.  
  784.   /* Try to account for error in evaluation of sin(y), cos(y).
  785.    * See above.
  786.    * FIXME ?
  787.    */
  788.   if(y > 1.0/GSL_DBL_EPSILON) {
  789.     result->err *= 0.5 * y;
  790.   }
  791.   else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
  792.     result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
  793.   }
  794.  
  795.   return GSL_SUCCESS;
  796. }
  797.  
  798.  
  799. /************************************************************************
  800.  *                                                                      *
  801.   Asymptotic approximations 8.11.5, 8.12.5, and 8.42.7 from
  802.   G.N.Watson, A Treatise on the Theory of Bessel Functions,
  803.   2nd Edition (Cambridge University Press, 1944).
  804.   Higher terms in expansion for x near l given by
  805.   Airey in Phil. Mag. 31, 520 (1916).
  806.  
  807.   This approximation is accurate to near 0.1% at the boundaries
  808.   between the asymptotic regions; well away from the boundaries
  809.   the accuracy is better than 10^{-5}.
  810.  *                                                                      *
  811.  ************************************************************************/
  812. #if 0
  813. double besselJ_meissel(double nu, double x)
  814. {
  815.   double beta = pow(nu, 0.325);
  816.   double result;
  817.  
  818.   /* Fitted matching points.   */
  819.   double llimit = 1.1 * beta;
  820.   double ulimit = 1.3 * beta;
  821.  
  822.   double nu2 = nu * nu;
  823.  
  824.   if (nu < 5. && x < 1.)
  825.     {
  826.       /* Small argument and order. Use a Taylor expansion. */
  827.       int k;
  828.       double xo2 = 0.5 * x;
  829.       double gamfactor = pow(nu,nu) * exp(-nu) * sqrt(nu * 2. * M_PI)
  830.     * (1. + 1./(12.*nu) + 1./(288.*nu*nu));
  831.       double prefactor = pow(xo2, nu) / gamfactor;
  832.       double C[5];
  833.  
  834.       C[0] = 1.;
  835.       C[1] = -C[0] / (nu+1.);
  836.       C[2] = -C[1] / (2.*(nu+2.));
  837.       C[3] = -C[2] / (3.*(nu+3.));
  838.       C[4] = -C[3] / (4.*(nu+4.));
  839.       
  840.       result = 0.;
  841.       for(k=0; k<5; k++)
  842.     result += C[k] * pow(xo2, 2.*k);
  843.  
  844.       result *= prefactor;
  845.     }
  846.   else if(x < nu - llimit)
  847.     {
  848.       /* Small x region: x << l.    */
  849.       double z = x / nu;
  850.       double z2 = z*z;
  851.       double rtomz2 = sqrt(1.-z2);
  852.       double omz2_2 = (1.-z2)*(1.-z2);
  853.  
  854.       /* Calculate Meissel exponent. */
  855.       double term1 = 1./(24.*nu) * ((2.+3.*z2)/((1.-z2)*rtomz2) -2.);
  856.       double term2 = - z2*(4. + z2)/(16.*nu2*(1.-z2)*omz2_2);
  857.       double V_nu = term1 + term2;
  858.       
  859.       /* Calculate the harmless prefactor. */
  860.       double sterlingsum = 1. + 1./(12.*nu) + 1./(288*nu2);
  861.       double harmless = 1. / (sqrt(rtomz2*2.*M_PI*nu) * sterlingsum);
  862.  
  863.       /* Calculate the logarithm of the nu dependent prefactor. */
  864.       double ln_nupre = rtomz2 + log(z) - log(1. + rtomz2);
  865.  
  866.       result = harmless * exp(nu*ln_nupre - V_nu);
  867.     } 
  868.   else if(x < nu + ulimit)
  869.     {         
  870.       /* Intermediate region 1: x near nu. */
  871.       double eps = 1.-nu/x;
  872.       double eps_x = eps * x;
  873.       double eps_x_2 = eps_x * eps_x;
  874.       double xo6 = x/6.;
  875.       double B[6];
  876.       static double gam[6] = {2.67894, 1.35412, 1., 0.89298, 0.902745, 1.};
  877.       static double sf[6] = {0.866025, 0.866025, 0., -0.866025, -0.866025, 0.};
  878.       
  879.       /* Some terms are identically zero, because sf[] can be zero.
  880.        * Some terms do not appear in the result.
  881.        */
  882.       B[0] = 1.;
  883.       B[1] = eps_x;
  884.       /* B[2] = 0.5 * eps_x_2 - 1./20.; */
  885.       B[3] = eps_x * (eps_x_2/6. - 1./15.);
  886.       B[4] = eps_x_2 * (eps_x_2 - 1.)/24. + 1./280.;
  887.       /* B[5] = eps_x * (eps_x_2*(0.5*eps_x_2 - 1.)/60. + 43./8400.); */
  888.  
  889.       result  = B[0] * gam[0] * sf[0] / pow(xo6, 1./3.);
  890.       result += B[1] * gam[1] * sf[1] / pow(xo6, 2./3.);
  891.       result += B[3] * gam[3] * sf[3] / pow(xo6, 4./3.);
  892.       result += B[4] * gam[4] * sf[4] / pow(xo6, 5./3.);
  893.  
  894.       result /= (3.*M_PI);
  895.     }
  896.   else 
  897.     {
  898.       /* Region of very large argument. Use expansion
  899.        * for x>>l, and we need not be very exacting.
  900.        */
  901.       double secb = x/nu;
  902.       double sec2b= secb*secb;
  903.       
  904.       double cotb = 1./sqrt(sec2b-1.);      /* cotb=cot(beta) */
  905.  
  906.       double beta = acos(nu/x);
  907.       double trigarg = nu/cotb - nu*beta - 0.25 * M_PI;
  908.       
  909.       double cot3b = cotb * cotb * cotb;
  910.       double cot6b = cot3b * cot3b;
  911.  
  912.       double sum1, sum2, expterm, prefactor, trigcos;
  913.  
  914.       sum1  = 2.0 + 3.0 * sec2b;
  915.       trigarg -= sum1 * cot3b / (24.0 * nu);
  916.  
  917.       trigcos = cos(trigarg);
  918.  
  919.       sum2 = 4.0 + sec2b;
  920.       expterm = sum2 * sec2b * cot6b / (16.0 * nu2);
  921.  
  922.       expterm = exp(-expterm);
  923.       prefactor = sqrt(2. * cotb / (nu * M_PI));
  924.       
  925.       result = prefactor * expterm * trigcos;
  926.     }
  927.  
  928.   return  result;
  929. }
  930. #endif
  931.